<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB R2018b"><title>Thick Plate Finite Element Model</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.7999992370605px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 13.9999990463257px; font-weight: normal; text-align: left;  }
.S2 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S3 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 13.9999990463257px; font-weight: normal; text-align: center;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S4 { border-left: 1.44000005722046px solid rgb(233, 233, 233); border-right: 1.44000005722046px solid rgb(233, 233, 233); border-top: 1.44000005722046px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 13.9999990463257px;  }
.S5 { border-left: 1.44000005722046px solid rgb(233, 233, 233); border-right: 1.44000005722046px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 13.9999990463257px;  }
.S6 { border-left: 1.44000005722046px solid rgb(233, 233, 233); border-right: 1.44000005722046px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1.44000005722046px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 13.9999990463257px;  }
.S7 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 13.9999990463257px; overflow-x: hidden; line-height: 17.234001159668px;  }
.S8 { border-left: 1.44000005722046px solid rgb(233, 233, 233); border-right: 1.44000005722046px solid rgb(233, 233, 233); border-top: 1.44000005722046px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 13.9999990463257px;  }
.embeddedOutputsErrorElement {min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
.embeddedOutputsWarningElement{min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
.diagnosticMessage-wrapper {font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {color: rgb(255,100,0); text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {color: rgb(230,0,0); text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart {white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {white-space: pre;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {white-space: pre; word-wrap: initial; min-height: 18px; max-height: 250px; overflow: auto;}
.textElement,.rtcDataTipElement .textElement {padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {min-height: 16px;}
.rightPaneElement .textElement {padding-top: 2px; padding-left: 9px;}
.S9 { border-left: 1.44000005722046px solid rgb(233, 233, 233); border-right: 1.44000005722046px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1.44000005722046px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 13.9999990463257px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Thick Plate Finite Element Model</span></h1><div  class = 'S1'><span>by Simone Coniglio and prof Joseph Morlier ISAE-SUPAERO</span></div><h2  class = 'S2'><span>Geometry interpolation</span></h2><div  class = 'S1'><span>Gemetric description and displacement interpolation of a Quad4 finite element</span></div><div  class = 'S1'><span>Given the node coordinates:</span></div><div  class = 'S3'><span style="vertical-align:-29px"><img src="" width="139" height="66" /></span></div><div  class = 'S1'><span>the middle plane geometry can be interpolated throught bilinear interpolation:</span></div><div  class = 'S3'><span style="vertical-align:-15px"><img src="" width="288.5" height="34.5" /></span></div><div  class = 'S1'><span>Where the 4 bilinear shape function are:</span></div><div  class = 'S3'><span style="vertical-align:-15px"><img src="" width="183.5" height="31.5" /></span></div><div  class = 'S1'><span>and </span><span style="vertical-align:-6px"><img src="" width="23.5" height="20" /></span><span> are the vertex normal.</span></div><div  class = 'S1'><span>To find this vector we can compute the vectors:</span></div><div  class = 'S3'><span style="vertical-align:-6px"><img src="" width="37" height="20" /></span><span style="vertical-align:-15px"><img src="" width="128.5" height="36" /></span><span> </span></div><div  class = 'S3'><span style="vertical-align:-6px"><img src="" width="37.5" height="20" /></span><span style="vertical-align:-15px"><img src="" width="131" height="36" /></span></div><div  class = 'S1'><span></span></div><div  class = 'S1'><span>Then  </span><span style="vertical-align:-6px"><img src="" width="25" height="23" /></span><span> is:</span></div><div  class = 'S3'><span style="vertical-align:-17px"><img src="" width="114.5" height="39" /></span></div><h2  class = 'S2'><span>Displacement interpolation and strain-displacement matrix</span></h2><div  class = 'S1'><span>The only hypothesys made for this kind of element is the fact that the thickness is rigid and the consequent stress vanishes.</span></div><div  class = 'S1'><span>The displacement will be considered as:</span></div><div  class = 'S3'><span style="vertical-align:-15px"><img src="" width="397" height="34.5" /></span></div><div  class = 'S1'><span>It is possible to define distinct shape functions for angles and displacements. If we consentrate on the transversal behavior:</span></div><div  class = 'S3'><span style="vertical-align:-5px"><img src="" width="27.5" height="17.5" /></span><span style="vertical-align:-11px"><img src="" width="105" height="30" /></span></div><div  class = 'S3'><span style="vertical-align:-5px"><img src="" width="27.5" height="17.5" /></span><span style="vertical-align:-11px"><img src="" width="108" height="30" /></span></div><div  class = 'S1'><span>where </span></div><div  class = 'S3'><span style="vertical-align:-17px"><img src="" width="133" height="40.5" /></span><span style="vertical-align:-17px"><img src="" width="30.5" height="40.5" /></span></div><div  class = 'S1'><span>Defining</span></div><div  class = 'S1'><span></span></div><div  class = 'S3'><span style="vertical-align:-5px"><img src="" width="17" height="17.5" /></span><span style="vertical-align:-51px"><img src="" width="76.5" height="112.5" /></span></div><div  class = 'S1'><span>so that we can call:</span></div><div  class = 'S3'><span style="vertical-align:-6px"><img src="" width="31.5" height="20" /></span><span style="vertical-align:-17px"><img src="" width="218.5" height="42" /></span></div><div  class = 'S3'><span style="vertical-align:-137px"><img src="" width="432.5" height="282" /></span></div><div  class = 'S3'><span style="vertical-align:-6px"><img src="" width="94" height="20" /></span><span style="vertical-align:-11px"><img src="" width="165.5" height="30" /></span><span style="vertical-align:-14px"><img src="" width="112" height="37.5" /></span></div><div  class = 'S3'><span style="vertical-align:-137px"><img src="" width="430.5" height="282" /></span></div><div  class = 'S3'><span style="vertical-align:-89px"><img src="" width="449.5" height="187" /></span></div><h2  class = 'S2'><span>Momentum/Shear - deformation relations</span></h2><div  class = 'S1'><span>Defining the material property:</span></div><div  class = 'S3'><span style="vertical-align:-16px"><img src="" width="96.5" height="33" /></span></div><div  class = 'S1'><span>then one can compute the moments vector as:</span></div><div  class = 'S3'><span style="vertical-align:-33px"><img src="" width="226" height="76.5" /></span></div><div  class = 'S1'><span>and the shear forces as:</span></div><div  class = 'S3'><span style="vertical-align:-17px"><img src="" width="154.5" height="40.5" /></span></div><div  class = 'S1'><span>for the membrane behavior:</span></div><div  class = 'S3'><span style="vertical-align:-29px"><img src="" width="119.5" height="66" /></span><span style="vertical-align:-33px"><img src="" width="142" height="76.5" /></span></div><h2  class = 'S2'><span>Local/global displacements</span></h2><div  class = 'S1'><span>First of all we will use the several convention for the local DOFs:</span></div><div  class = 'S3'><span style="vertical-align:-6px"><img src="" width="134.5" height="20" /></span></div><div  class = 'S1'><span>where</span></div><div  class = 'S3'><span style="vertical-align:-6px"><img src="" width="169" height="20" /></span></div><div  class = 'S1'><span>the Global DOF vector can be expressed as:</span></div><div  class = 'S3'><span style="vertical-align:-6px"><img src="" width="190.5" height="26" /></span></div><div  class = 'S3'><span style="vertical-align:-6px"><img src="" width="36.5" height="20" /></span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: normal; font-weight: normal;">=</span><span style="vertical-align:-10px"><img src="" width="72.5" height="34.5" /></span></div><div  class = 'S1'><span></span></div><h2  class = 'S2'><span>Stiffness Matrix</span></h2><div  class = 'S3'><span style="vertical-align:-12px"><img src="" width="199.5" height="33" /></span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre;"><span>nx=10;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>ny=10;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>nu=0.3;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>t=1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>E=1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>R=10;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>theta=linspace(0,pi/2,nx+1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>zeta=linspace(0,10,ny+1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>[Tet,Z]=meshgrid(theta,zeta);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>COORD=[R*cos(Tet(:)');</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    R*sin(Tet(:)');</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    Z(:)']';</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>nodenrs=reshape(1:((nx+1)*(ny+1)),ny+1,nx+1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>ELEMENT=zeros(nx*ny,4);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>elid=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">for </span><span>k=1:ny</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span style="color: rgb(0, 0, 255);">for </span><span>l=1:nx</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        elid=elid+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        ELEMENT(elid,:)=[nodenrs(k,l) nodenrs(k+1,l) nodenrs(k+1,l+1) nodenrs(k,l+1)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span style="color: rgb(0, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>ELEMENT_DOF=[6*(ELEMENT(:,1)-1)+1,6*(ELEMENT(:,1)-1)+2,6*(ELEMENT(:,1)-1)+3,6*(ELEMENT(:,1)-1)+4,6*(ELEMENT(:,1)-1)+5,6*(ELEMENT(:,1)-1)+6];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>ELEMENT_DOF=[ELEMENT_DOF,6*(ELEMENT(:,2)-1)+1,6*(ELEMENT(:,2)-1)+2,6*(ELEMENT(:,2)-1)+3,6*(ELEMENT(:,2)-1)+4,6*(ELEMENT(:,2)-1)+5,6*(ELEMENT(:,2)-1)+6];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>ELEMENT_DOF=[ELEMENT_DOF,6*(ELEMENT(:,3)-1)+1,6*(ELEMENT(:,3)-1)+2,6*(ELEMENT(:,3)-1)+3,6*(ELEMENT(:,3)-1)+4,6*(ELEMENT(:,3)-1)+5,6*(ELEMENT(:,3)-1)+6];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>ELEMENT_DOF=[ELEMENT_DOF,6*(ELEMENT(:,4)-1)+1,6*(ELEMENT(:,4)-1)+2,6*(ELEMENT(:,4)-1)+3,6*(ELEMENT(:,4)-1)+4,6*(ELEMENT(:,4)-1)+5,6*(ELEMENT(:,4)-1)+6];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>figure</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>patch(</span><span style="color: rgb(160, 32, 240);">'Faces'</span><span>,ELEMENT,</span><span style="color: rgb(160, 32, 240);">'Vertices'</span><span>,COORD,</span><span style="color: rgb(160, 32, 240);">'FaceColor'</span><span>,</span><span style="color: rgb(160, 32, 240);">'c'</span><span>)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>axis </span><span style="color: rgb(160, 32, 240);">equal</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: pre;"><span>view(-45,24)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="78D72921" data-scroll-top="null" data-scroll-left="null" data-testid="output_0" style="width: 1250px;"><div class="figureElement"><img class="figureImage figureContainingNode" src=""></div></div></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre;"><span>gauss_point=1/sqrt(3)*[-1 1];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>[csi,eta]=meshgrid(gauss_point,gauss_point);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N1=@(a,b) 1/4*(1-a)*(1-b);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N2=@(a,b) 1/4*(1+a)*(1-b);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N3=@(a,b) 1/4*(1+a)*(1+b);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N4=@(a,b) 1/4*(1-a)*(1+b);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N1_a=@(a,b) -1/4*(1-b);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N1_b=@(a,b) -1/4*(1-a);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N2_a=@(a,b) 1/4*(1-b);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N2_b=@(a,b) -1/4*(1+a);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N3_a=@(a,b) 1/4*(1+b);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N3_b=@(a,b) 1/4*(1+a);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N4_a=@(a,b) -1/4*(1+b);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>N4_b=@(a,b) 1/4*(1-a);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>B=zeros(8,24,4);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">% nodal locations</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>x1=COORD(ELEMENT(:,1),:);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>x2=COORD(ELEMENT(:,2),:);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>x3=COORD(ELEMENT(:,3),:);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>x4=COORD(ELEMENT(:,4),:);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">% element coordinate system</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>x0=1/4*(x1+x2+x3+x4);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>xE=1/2*(x2+x3)-x0; xE=xE./repmat(sqrt(sum(xE.^2,2)),1,3);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>zE=cross3(x2-x0,x3-x0); zE=zE./repmat(sqrt(sum(zE.^2,2)),1,3);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>yE = cross3(zE,xE);yE=yE./repmat(sqrt(sum(yE.^2,2)),1,3);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">% loop over element for stiffness matrix assembly</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>I=[];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>J=[];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>Kij=[];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>Eid=[];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">for </span><span>k=1:size(x0,1)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    TEG = [xE(k,:)',yE(k,:)',zE(k,:)'].';</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span style="color: rgb(34, 139, 34);">% node positions in element coordinate system</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    XE = TEG *([x1(k,:)',x2(k,:)',x3(k,:)',x4(k,:)']-[x0(k,:)',x0(k,:)',x0(k,:)',x0(k,:)']);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    obj.XE=XE;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span style="color: rgb(34, 139, 34);">% nodal unit normals</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    j1=@(xi,eta) .25*XE*[      -(1-eta),       (1-eta),       (1+eta),      -(1+eta)].';</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    j2=@(xi,eta) .25*XE*[(1-xi)*-1     ,(1+xi)*-1     ,(1+xi)        ,(1-xi)        ].';</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    obj.n1 = cross3(j1(-1,-1)',j2(-1,-1)'); obj.n1 = obj.n1./norm(obj.n1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    obj.n2 = cross3(j1( 1,-1)',j2( 1,-1)'); obj.n2 = obj.n2./norm(obj.n2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    obj.n3 = cross3(j1( 1, 1)',j2( 1, 1)'); obj.n3 = obj.n3./norm(obj.n3);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    obj.n4 = cross3(j1(-1, 1)',j2(-1, 1)'); obj.n4 = obj.n4./norm(obj.n4);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    obj.G= [t/(1-nu^2)*[1 nu 0;nu 1 0;0 0 (1-nu)/2]          ,zeros(3,3)                     , zeros(3,2) </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        zeros(3,3)  ,t^3/12/(1-nu^2)*[1 nu 0;nu 1 0;0 0 (1-nu)/2]  , zeros(3,2)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        zeros(2,3)  ,zeros(2,3)                     ,5/6*t/2/(1+nu)*eye(2)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">    %% Coordinate transformation at nodes</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>obj.A = zeros(3,3,4);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">for </span><span>n =1:4</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span style="color: rgb(0, 0, 255);">switch </span><span>n</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span><span style="color: rgb(0, 0, 255);">case </span><span>1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            z_n = obj.n1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span><span style="color: rgb(0, 0, 255);">case </span><span>2</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            z_n = obj.n2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span><span style="color: rgb(0, 0, 255);">case </span><span>3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            z_n = obj.n3;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span><span style="color: rgb(0, 0, 255);">case </span><span>4</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            z_n = obj.n4;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span style="color: rgb(0, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    x_n = cross3([0 1 0],z_n)./norm(cross3([0 1 0],z_n));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    y_n = cross3(z_n,x_n);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    TNIE= [x_n;y_n;z_n];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    obj.A(1:3,1:3,n) = TNIE.'*[0 1 0; -1 0 0; 0 0 1]*TNIE;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>XI = 1/sqrt(3)*[-1 -1; 1 1];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>ETA = 1/sqrt(3)*[-1 1; -1 1];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>ke=zeros(24);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">for </span><span>i = 1:2</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span style="color: rgb(0, 0, 255);">for </span><span>j = 1:2</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        xi = XI(i,j);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        eta = ETA(i,j);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">        %% shape function evaluations</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        Ni      = .25*[(1-xi)*(1-eta),(1+xi)*(1-eta),(1+xi)*(1+eta),(1-xi)*(1+eta)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        dNdxii  = .25*[      -(1-eta),       (1-eta),       (1+eta),      -(1+eta)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        dNdetai = .25*[(1-xi)*-1     ,(1+xi)*-1     ,(1+xi)        ,(1-xi)        ];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">        %% Jacobian</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        Jac = </span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            [ [dNdxii;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdetai]* obj.XE.';</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            Ni*[obj.n1',obj.n2',obj.n3',obj.n4'].' ];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        detJ = det(Jac);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">        %% Rotation matrix</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z_i = cross3(Jac(1,:),Jac(2,:))./norm(cross3(Jac(1,:),Jac(2,:)));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        x_i = cross3([0 1 0],z_i)./norm(cross3([0 1 0],z_i));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        y_i = cross3(z_i,x_i);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        T = [x_i;y_i;z_i];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">        %% Partial derivatives wrt physical coordinates</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        dNdxzy = T/Jac*[dNdxii; dNdetai; 0 0 0 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">        %% Element thickness at integration point</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        tg = t;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span><span style="color: rgb(34, 139, 34);">% Calculate strain displacement matrix at xi &amp; eta</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        e1=[dNdxzy(1,1)     0           0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0               dNdxzy(2,1) 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdxzy(2,1)     dNdxzy(1,1) 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        e2=[dNdxzy(1,2)     0           0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0               dNdxzy(2,2) 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdxzy(2,2)     dNdxzy(1,2) 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        e3=[dNdxzy(1,3)     0           0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0               dNdxzy(2,3) 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdxzy(2,3)     dNdxzy(1,3) 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        e4=[dNdxzy(1,4)     0           0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0               dNdxzy(2,4) 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdxzy(2,4)     dNdxzy(1,4) 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s1=[0 0 dNdxzy(2,1) 0    Ni(1) 0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0 0 dNdxzy(1,1) Ni(1) 0    0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s2=[0 0 dNdxzy(2,2) 0    Ni(2) 0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0 0 dNdxzy(1,2) Ni(2) 0    0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s3=[0 0 dNdxzy(2,3) 0    Ni(3) 0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0 0 dNdxzy(1,3) Ni(3) 0    0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s4=[0 0 dNdxzy(2,4) 0    Ni(4) 0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0 0 dNdxzy(1,4) Ni(4) 0    0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3 = zeros(3);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>         B  = [[e1  z3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3  e1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s1                  ] * [T z3; z3 T*obj.A(:,:,1)],</span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        [e2  z3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3  e2</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s2                  ] * [T z3; z3 T*obj.A(:,:,2)],</span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        [e3  z3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3  e3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s3                  ] * [T z3; z3 T*obj.A(:,:,3)],</span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        [e4  z3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3  e4</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s4                  ] * [T z3; z3 T*obj.A(:,:,4)]];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">    %% Modify constitutive (stress-strain) matrix for bending</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    G = obj.G;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    ke =ke+ B([1:2,4:5],:).'*G([1:2,4:5],[1:2,4:5])*B([1:2,4:5],:) *tg*detJ;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span style="color: rgb(0, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span> </span><span style="color: rgb(34, 139, 34);">% 1 point integration for shear</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span> xi=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span> eta=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">  %% shape function evaluations</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        Ni      = .25*[(1-xi)*(1-eta),(1+xi)*(1-eta),(1+xi)*(1+eta),(1-xi)*(1+eta)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        dNdxii  = .25*[      -(1-eta),       (1-eta),       (1+eta),      -(1+eta)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        dNdetai = .25*[(1-xi)*-1     ,(1+xi)*-1     ,(1+xi)        ,(1-xi)        ];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">        %% Jacobian</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        Jac = </span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            [ [dNdxii;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdetai]* obj.XE.';</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            Ni*[obj.n1',obj.n2',obj.n3',obj.n4'].' ];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        detJ = det(Jac);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">        %% Rotation matrix</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z_i = cross3(Jac(1,:),Jac(2,:))./norm(cross3(Jac(1,:),Jac(2,:)));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        x_i = cross3([0 1 0],z_i)./norm(cross3([0 1 0],z_i));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        y_i = cross3(z_i,x_i);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        T = [x_i;y_i;z_i];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">        %% Partial derivatives wrt physical coordinates</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        dNdxzy = T/Jac*[dNdxii; dNdetai; 0 0 0 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">        %% Element thickness at integration point</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        tg = t;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span><span style="color: rgb(34, 139, 34);">% Calculate strain displacement matrix at xi &amp; eta</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        e1=[dNdxzy(1,1)     0           0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0               dNdxzy(2,1) 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdxzy(2,1)     dNdxzy(1,1) 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        e2=[dNdxzy(1,2)     0           0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0               dNdxzy(2,2) 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdxzy(2,2)     dNdxzy(1,2) 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        e3=[dNdxzy(1,3)     0           0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0               dNdxzy(2,3) 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdxzy(2,3)     dNdxzy(1,3) 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        e4=[dNdxzy(1,4)     0           0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0               dNdxzy(2,4) 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            dNdxzy(2,4)     dNdxzy(1,4) 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s1=[0 0 dNdxzy(2,1) 0    Ni(1) 0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0 0 dNdxzy(1,1) Ni(1) 0    0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s2=[0 0 dNdxzy(2,2) 0    Ni(2) 0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0 0 dNdxzy(1,2) Ni(2) 0    0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s3=[0 0 dNdxzy(2,3) 0    Ni(3) 0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0 0 dNdxzy(1,3) Ni(3) 0    0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s4=[0 0 dNdxzy(2,4) 0    Ni(4) 0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>            0 0 dNdxzy(1,4) Ni(4) 0    0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3 = zeros(3);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>         B  = [[e1  z3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3  e1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s1                  ] * [T z3; z3 T*obj.A(:,:,1)],</span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        [e2  z3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3  e2</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s2                  ] * [T z3; z3 T*obj.A(:,:,2)],</span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        [e3  z3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3  e3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s3                  ] * [T z3; z3 T*obj.A(:,:,3)],</span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        [e4  z3</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        z3  e4</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>        s4                  ] * [T z3; z3 T*obj.A(:,:,4)]];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">    %% Modify constitutive (stress-strain) matrix for bending</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    G = obj.G;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    ks = B([3,6:8]  ,:).'*G([3,6:8]  ,[3,6:8]  )*B([3,6:8]  ,:) *tg*detJ;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    ke = ke + 4*ks;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">    %% transform to global coordinate system and save</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    REG(22:24,22:24)=TEG;REG(19:21,19:21)=TEG;REG(16:18,16:18)=TEG;REG(13:15,13:15)=TEG;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    REG(10:12,10:12)=TEG;REG(7:9,7:9)=TEG;REG(4:6,4:6)=TEG;REG(1:3,1:3)=TEG;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    RGE = REG.';</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span style="color: rgb(34, 139, 34);">% save the I,J,K for assembly</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    [ie,je,kke]=find(RGE*ke*REG);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    eid=k*ones(size(ie));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span class="warning_squiggle_rte">I</span><span>=[I;ELEMENT_DOF(k,ie)'];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span class="warning_squiggle_rte">J</span><span>=[J;ELEMENT_DOF(k,je)'];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span class="warning_squiggle_rte">Kij</span><span>=[Kij;kke];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    </span><span class="warning_squiggle_rte">Eid</span><span>=[Eid;eid];</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">end</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="E3F190A6" data-scroll-top="null" data-scroll-left="null" data-testid="output_1" style="width: 1250px;"><div class="figureElement"><img class="figureImage figureContainingNode" src=""></div></div></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre;"><span>clamped_nodes=find(COORD(:,3)==0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>clamped_DOF=[6*(clamped_nodes-1)+1;6*(clamped_nodes-1)+2;6*(clamped_nodes-1)+3;6*(clamped_nodes-1)+4;6*(clamped_nodes-1)+5;6*(clamped_nodes-1)+6];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>load_nodes=find(COORD(:,3)==10);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>extr_load_nodes=find(COORD(:,3)==10&amp;(COORD(:,2)==10|COORD(:,1)==10));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>extr_load_dof=</span><span class="warning_squiggle_rte">[</span><span>6*(extr_load_nodes-1)+3];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>load_DOF=</span><span class="warning_squiggle_rte">[</span><span>6*(load_nodes-1)+3];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>all_DOFs=(1:6*size(COORD,1))';</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>free_DOFs=setdiff(all_DOFs,clamped_DOF);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>load_module=[ones(size(load_DOF));-0.5*[1;1]];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>F=sparse([load_DOF;extr_load_dof],1,load_module,6*size(COORD,1),1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>K=sparse(I,J,E*Kij,6*size(COORD,1),6*size(COORD,1));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>K=(K+K')/2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>U=zeros(size(F));</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: pre;"><span>U(free_DOFs)=K(free_DOFs,free_DOFs)\F(free_DOFs);</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsWarningElement" uid="8F8D8675" data-scroll-top="null" data-scroll-left="null" data-width="1220" data-height="18" data-hashorizontaloverflow="false" data-testid="output_2" style="max-height: 261px; width: 1250px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is singular to working precision.</div><div class="diagnosticMessage-stackPart"></div></div></div></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre;"><span>figure</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>patch(</span><span style="color: rgb(160, 32, 240);">'Faces'</span><span>,ELEMENT,</span><span style="color: rgb(160, 32, 240);">'Vertices'</span><span>,COORD,</span><span style="color: rgb(160, 32, 240);">'FaceVertexCData'</span><span>,U(3:6:end),</span><span style="color: rgb(160, 32, 240);">'FaceColor'</span><span>,</span><span style="color: rgb(160, 32, 240);">'interp'</span><span>);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>colorbar</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>axis </span><span style="color: rgb(160, 32, 240);">equal</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre;"><span>view(-45,24)</span></span></div></div></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">function </span><span>p = cross3(u,v)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre;"><span>    p = [u(:,2).*v(:,3) u(:,3).*v(:,1) u(:,1).*v(:,2)]-[u(:,3).*v(:,2) u(:,1).*v(:,3) u(:,2).*v(:,1)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">end</span></span></div></div></div></div><br>
<!-- 
##### SOURCE BEGIN #####
%% Thick Plate Finite Element Model
% by Simone Coniglio and prof Joseph Morlier ISAE-SUPAERO
%% Geometry interpolation
% Gemetric description and displacement interpolation of a Quad4 finite element
% 
% Given the node coordinates:
% 
% $${\mathbf{X}}_{\mathbf{c}} =\left\lbrack \begin{array}{c}x_1  & x_2  & 
% x_3  & x_4 \\y_1  & y_2  & y_3  & y_4 \\z_1  & z_2  & z_3  & z_4 \end{array}\right\rbrack$$
% 
% the middle plane geometry can be interpolated throught bilinear interpolation:
% 
% $$\mathbf{X}\left(s,r,t\right)=\sum_{h=1}^{4} {\mathbf{X}}_{\mathbf{ch}} 
% N_h \left(s,r\right)+\frac{t}{2}a_h N_h \left(s,r\right){\hat{\mathbf{V}} }_{\mathbf{nh}}$$
% 
% Where the 4 bilinear shape function are:
% 
% $$N_h \left(s,r\right)=\frac{1}{4}\left(1+s_h s\right)\left(1+r_h r\right)$$
% 
% and ${\mathbf{V}}_{n\mathbf{\text{h}}}$ are the vertex normal.
% 
% To find this vector we can compute the vectors:
% 
% $${\mathbf{V}}_{s\mathbf{\text{h}}} =$$$$\sum_{h=1}^{4} {\mathbf{X}}_{\mathbf{ch}} 
% \frac{\partial N_h }{\partial s}\left(s_h ,r_h \right)$$ 
% 
% $${\mathbf{V}}_{\mathbf{rh}} =$$$$\sum_{h=1}^{4} {\mathbf{X}}_{\mathbf{\text{ch}}} 
% \frac{\partial N_h }{\partial r}\left(s_h ,r_h \right)$$
% 
% 
% 
% Then  ${\hat{\mathbf{V}} }_{n\mathbf{h}}$ is:
% 
% $${\hat{\mathbf{V}} }_{n\mathbf{h}} =\frac{{\mathbf{V}}_{s\mathbf{\text{h}}} 
% \times {\mathbf{V}}_{\mathbf{rh}} }{\left\|{\mathbf{V}}_{s\mathbf{\text{h}}} 
% \times {\mathbf{V}}_{\mathbf{rh}} \right\|}$$
%% Displacement interpolation and strain-displacement matrix
% The only hypothesys made for this kind of element is the fact that the thickness 
% is rigid and the consequent stress vanishes.
% 
% The displacement will be considered as:
% 
% $$\mathbf{u}\left(s,r,t\right)=\sum_{h=1}^{4} {\mathbf{u}}_{\mathbf{h}} 
% N_h \left(s,r\right)+\frac{t}{2}a_h N_{\theta h} \left(s,r\right)\theta_h {\times 
% \hat{\mathbf{V}} }_{\mathbf{nh}} =\left\lbrack H\left(s,r,t\right)\right\rbrack 
% \mathbf{q}$$
% 
% It is possible to define distinct shape functions for angles and displacements. 
% If we consentrate on the transversal behavior:
% 
% $$w=$$$$\sum_{h=1}^{4} w_{\mathbf{h}} N_h \left(s,r\right)$$
% 
% $$\phi =$$$$\sum_{h=1}^{4} {\phi \text{\,}}_{\mathbf{h}} N_h \left(s,r\right)$$
% 
% where 
% 
% $$\phi =\left\lbrack \begin{array}{c}\phi_s \\\phi_r \end{array}\right\rbrack 
% =\left\lbrack \begin{array}{c}0 & 1\\-1 & 0\end{array}\right\rbrack$$$$\left\lbrack 
% \begin{array}{c}\theta_s \\\theta_r \end{array}\right\rbrack$$
% 
% Defining
% 
% 
% 
% $$\mathcal{L}$$$$=\left\lbrack \begin{array}{c}\frac{\partial }{\partial 
% s} & 0\\0 & \frac{\partial }{\partial r}\\\frac{\partial }{\partial r} & \frac{\partial 
% }{\partial s}\end{array}\right\rbrack$$
% 
% so that we can call:
% 
% $$\varepsilon_m =$$$$\mathcal{L}\phi=\sum_{h=1}^4\mathcal{L}N_h \phi_h=\begin{array}{cc}[0 
% & B_{m\phi}]\end{array}\left\lbrace\begin{array}{c}w \\ \phi\end{array}\right\rbrace$$
% 
% $$\varepsilon_m =\left\lbrack \begin{array}{c}0 & 0 & 0 & 0 & \frac{\partial 
% N_1 }{\partial s} & 0 & \frac{\partial N_2 }{\partial s} & 0 & \frac{\partial 
% N_3 }{\partial s} & 0 & \frac{\partial N_4 }{\partial s} & 0\\0 & 0 & 0 & 0 
% & 0 & \frac{\partial N_1 }{\partial r} & 0 & \frac{\partial N_2 }{\partial r} 
% & 0 & \frac{\partial N_3 }{\partial r} & 0 & \frac{\partial N_4 }{\partial r}\\0 
% & 0 & 0 & 0 & \frac{\partial N_1 }{\partial r} & \frac{\partial N_1 }{\partial 
% s} & \frac{\partial N_2 }{\partial r} & \frac{\partial N_2 }{\partial s} & \frac{\partial 
% N_3 }{\partial r} & \frac{\partial N_3 }{\partial s} & \frac{\partial N_4 }{\partial 
% r} & \frac{\partial N_4 }{\partial s}\end{array}\right\rbrack \left\lbrack \begin{array}{c}w_1 
% \\w_2 \\w_3 \\w_4 \\\phi_{s1} \\\phi_{r1} \\\phi_{s2} \\\phi_{r2} \\\phi_{s3} 
% \\\phi_{r3} \\\phi_{s4} \\\phi_{r4} \end{array}\right\rbrack$$
% 
% $$\varepsilon_t =\nabla w+\phi =$$$$\sum_{h=1}^{4} w_{\mathbf{h}} {\nabla 
% \text{\,}N}_h \left(s,r\right)+N_h \phi_h$$$$=\begin{array}{cc}[B_{sw} & B_{s\phi}]\end{array}\left\lbrace\begin{array}{c}w 
% \\ \phi\end{array}\right\rbrace$$
% 
% $$\varepsilon_t =\left\lbrack \begin{array}{c}\frac{\partial N_1 }{\partial 
% s} & \frac{\partial N_2 }{\partial s} & \frac{\partial N_3 }{\partial s} & \frac{\partial 
% N_4 }{\partial s} & N_1  & 0 & N_2  & 0 & N_3  & 0 & N_4  & 0\\\frac{\partial 
% N_1 }{\partial r} & \frac{\partial N_2 }{\partial r} & \frac{\partial N_3 }{\partial 
% r} & \frac{\partial N_4 }{\partial r} & 0 & N_1  & 0 & N_2  & 0 & N_3  & 0 & 
% N_4 \end{array}\right\rbrack \left\lbrack \begin{array}{c}w_1 \\w_2 \\w_3 \\w_4 
% \\\phi_{s1} \\\phi_{\text{r1}} \\\phi_{\text{s2}} \\\phi_{\text{r2}} \\\phi_{\text{s3}} 
% \\\phi_{\text{r3}} \\\phi_{\text{s4}} \\\phi_{\text{r4}} \end{array}\right\rbrack$$
% 
% $$\left\lbrack \begin{array}{c}\varepsilon_s \\\varepsilon_r \\\varepsilon_{\text{rs}} 
% \end{array}\right\rbrack =\left\lbrack B_{uv} \right\rbrack \left\lbrace q_{uv} 
% \right\rbrace =\left\lbrack \begin{array}{c}\frac{\partial N_1 }{\partial s} 
% & 0 & \frac{\partial N_2 }{\partial s} & 0 & \frac{\partial N_3 }{\partial s} 
% & 0 & \frac{\partial N_4 }{\partial s} & 0\\0 & \frac{\partial N_1 }{\partial 
% r} & 0 & \frac{\partial N_2 }{\partial r} & 0 & \frac{\partial N_3 }{\partial 
% r} & 0 & \frac{\partial N_4 }{\partial r}\\\frac{\partial N_1 }{\partial r} 
% & \frac{\partial N_1 }{\partial s} & \frac{\partial N_2 }{\partial r} & \frac{\partial 
% N_2 }{\partial s} & \frac{\partial N_3 }{\partial r} & \frac{\partial N_3 }{\partial 
% s} & \frac{\partial N_4 }{\partial r} & \frac{\partial N_4 }{\partial s}\end{array}\right\rbrack 
% \left\lbrack \begin{array}{c}u_1 \\v_1 \\u_2 \\v_2 \\u_3 \\v_3 \\u_4 \\v_4 \end{array}\right\rbrack$$
%% Momentum/Shear - deformation relations
% Defining the material property:
% 
% $$D=\frac{{Et}^3 }{12\left(1-\nu^2 \right)}$$
% 
% then one can compute the moments vector as:
% 
% $$\mathbf{M}=\left\lbrack \begin{array}{c}M_s \\M_r \\M_{sr} \end{array}\right\rbrack 
% =D\left\lbrack \begin{array}{c}1 & \nu  & 0\\\nu  & 1 & 0\\0 & 0 & \frac{\left(1-\nu 
% \text{\,}\right)}{2}\end{array}\right\rbrack \varepsilon_m$$
% 
% and the shear forces as:
% 
% $$\mathbf{T}=\left\lbrack \begin{array}{c}T_s \\T_r \end{array}\right\rbrack 
% =\kappa \text{Gt}\left\lbrack \begin{array}{c}1 & 0\\0 & 1\end{array}\right\rbrack 
% \varepsilon_t$$
% 
% for the membrane behavior:
% 
% $$\sigma =\left\lbrack \begin{array}{c}\sigma_s \\\sigma_r \\\sigma_{rs} 
% \end{array}\right\rbrack =\frac{E}{1-\nu^2 }$$$$\left\lbrack \begin{array}{c}1 
% & \nu  & 0\\\nu  & 1 & 0\\0 & 0 & \frac{\left(1-\nu \text{\,}\right)}{2}\end{array}\right\rbrack 
% \left\lbrack \begin{array}{c}\varepsilon_s \\\varepsilon_r \\\varepsilon_{rs} 
% \end{array}\right\rbrack$$
%% Local/global displacements
% First of all we will use the several convention for the local DOFs:
% 
% $$\left\lbrace q\right\rbrace ={\left\lbrace q_1 ,q_2 ,q_3 ,q_4 \right\rbrace 
% }^T$$
% 
% where
% 
% $$q_e =\left\lbrace u_e ,v_e ,w_e ,\phi_{se} ,\phi_{re} ,\phi_{te} \right\rbrace$$
% 
% the Global DOF vector can be expressed as:
% 
% $${\tilde{q} }_e =\left\lbrace u_{xe} ,v_{ye} ,w_{ze} ,\phi_{x\text{e}} 
% ,\phi_{y\text{e}} ,\phi_{z\text{e}} \right\rbrace$$
% 
% $${\left\lbrace q_e \right\rbrace }^T$$$$=$$$${\left\lbrack T_e \right\rbrack 
% }_{\text{l}g} {\left\lbrace \tilde{q_e } \right\rbrace }^T$$
% 
% 
%% Stiffness Matrix
% $$\left\lbrack K\right\rbrack =\int {{\left\lbrack T\right\rbrack }^T \left\lbrack 
% B\right\rbrack }^T \left\lbrack D\right\rbrack \left\lbrack B\right\rbrack \left\lbrack 
% T\right\rbrack d\Omega \text{\,}$$

nx=10;
ny=10;
nu=0.3;
t=1;
E=1;
R=10;
theta=linspace(0,pi/2,nx+1);
zeta=linspace(0,10,ny+1);
[Tet,Z]=meshgrid(theta,zeta);
COORD=[R*cos(Tet(:)');
    R*sin(Tet(:)');
    Z(:)']';
nodenrs=reshape(1:((nx+1)*(ny+1)),ny+1,nx+1);
ELEMENT=zeros(nx*ny,4);
elid=0;
for k=1:ny
    for l=1:nx
        elid=elid+1;
        ELEMENT(elid,:)=[nodenrs(k,l) nodenrs(k+1,l) nodenrs(k+1,l+1) nodenrs(k,l+1)];
    end
end
ELEMENT_DOF=[6*(ELEMENT(:,1)-1)+1,6*(ELEMENT(:,1)-1)+2,6*(ELEMENT(:,1)-1)+3,6*(ELEMENT(:,1)-1)+4,6*(ELEMENT(:,1)-1)+5,6*(ELEMENT(:,1)-1)+6];
ELEMENT_DOF=[ELEMENT_DOF,6*(ELEMENT(:,2)-1)+1,6*(ELEMENT(:,2)-1)+2,6*(ELEMENT(:,2)-1)+3,6*(ELEMENT(:,2)-1)+4,6*(ELEMENT(:,2)-1)+5,6*(ELEMENT(:,2)-1)+6];
ELEMENT_DOF=[ELEMENT_DOF,6*(ELEMENT(:,3)-1)+1,6*(ELEMENT(:,3)-1)+2,6*(ELEMENT(:,3)-1)+3,6*(ELEMENT(:,3)-1)+4,6*(ELEMENT(:,3)-1)+5,6*(ELEMENT(:,3)-1)+6];
ELEMENT_DOF=[ELEMENT_DOF,6*(ELEMENT(:,4)-1)+1,6*(ELEMENT(:,4)-1)+2,6*(ELEMENT(:,4)-1)+3,6*(ELEMENT(:,4)-1)+4,6*(ELEMENT(:,4)-1)+5,6*(ELEMENT(:,4)-1)+6];
figure
patch('Faces',ELEMENT,'Vertices',COORD,'FaceColor','c')
axis equal
view(-45,24)
gauss_point=1/sqrt(3)*[-1 1];
[csi,eta]=meshgrid(gauss_point,gauss_point);
N1=@(a,b) 1/4*(1-a)*(1-b);
N2=@(a,b) 1/4*(1+a)*(1-b);
N3=@(a,b) 1/4*(1+a)*(1+b);
N4=@(a,b) 1/4*(1-a)*(1+b);
N1_a=@(a,b) -1/4*(1-b);
N1_b=@(a,b) -1/4*(1-a);
N2_a=@(a,b) 1/4*(1-b);
N2_b=@(a,b) -1/4*(1+a);
N3_a=@(a,b) 1/4*(1+b);
N3_b=@(a,b) 1/4*(1+a);
N4_a=@(a,b) -1/4*(1+b);
N4_b=@(a,b) 1/4*(1-a);
B=zeros(8,24,4);
% nodal locations
x1=COORD(ELEMENT(:,1),:);
x2=COORD(ELEMENT(:,2),:);
x3=COORD(ELEMENT(:,3),:);
x4=COORD(ELEMENT(:,4),:);
% element coordinate system
x0=1/4*(x1+x2+x3+x4);
xE=1/2*(x2+x3)-x0; xE=xE./repmat(sqrt(sum(xE.^2,2)),1,3);
zE=cross3(x2-x0,x3-x0); zE=zE./repmat(sqrt(sum(zE.^2,2)),1,3);
yE = cross3(zE,xE);yE=yE./repmat(sqrt(sum(yE.^2,2)),1,3);
% loop over element for stiffness matrix assembly
I=[];
J=[];
Kij=[];
Eid=[];
for k=1:size(x0,1)
    TEG = [xE(k,:)',yE(k,:)',zE(k,:)'].';
    % node positions in element coordinate system
    XE = TEG *([x1(k,:)',x2(k,:)',x3(k,:)',x4(k,:)']-[x0(k,:)',x0(k,:)',x0(k,:)',x0(k,:)']);
    obj.XE=XE;
    % nodal unit normals
    j1=@(xi,eta) .25*XE*[      -(1-eta),       (1-eta),       (1+eta),      -(1+eta)].';
    j2=@(xi,eta) .25*XE*[(1-xi)*-1     ,(1+xi)*-1     ,(1+xi)        ,(1-xi)        ].';
    obj.n1 = cross3(j1(-1,-1)',j2(-1,-1)'); obj.n1 = obj.n1./norm(obj.n1);
    obj.n2 = cross3(j1( 1,-1)',j2( 1,-1)'); obj.n2 = obj.n2./norm(obj.n2);
    obj.n3 = cross3(j1( 1, 1)',j2( 1, 1)'); obj.n3 = obj.n3./norm(obj.n3);
    obj.n4 = cross3(j1(-1, 1)',j2(-1, 1)'); obj.n4 = obj.n4./norm(obj.n4);
    obj.G= [t/(1-nu^2)*[1 nu 0;nu 1 0;0 0 (1-nu)/2]          ,zeros(3,3)                     , zeros(3,2) 
        zeros(3,3)  ,t^3/12/(1-nu^2)*[1 nu 0;nu 1 0;0 0 (1-nu)/2]  , zeros(3,2)
        zeros(2,3)  ,zeros(2,3)                     ,5/6*t/2/(1+nu)*eye(2)];
    %% Coordinate transformation at nodes
obj.A = zeros(3,3,4);
for n =1:4
    switch n
        case 1
            z_n = obj.n1;
        case 2
            z_n = obj.n2;
        case 3
            z_n = obj.n3;
        case 4
            z_n = obj.n4;
    end
    x_n = cross3([0 1 0],z_n)./norm(cross3([0 1 0],z_n));
    y_n = cross3(z_n,x_n);
    TNIE= [x_n;y_n;z_n];
    obj.A(1:3,1:3,n) = TNIE.'*[0 1 0; -1 0 0; 0 0 1]*TNIE;
end
XI = 1/sqrt(3)*[-1 -1; 1 1];
ETA = 1/sqrt(3)*[-1 1; -1 1];
ke=zeros(24);
for i = 1:2
    for j = 1:2
        xi = XI(i,j);
        eta = ETA(i,j);
        %% shape function evaluations
        Ni      = .25*[(1-xi)*(1-eta),(1+xi)*(1-eta),(1+xi)*(1+eta),(1-xi)*(1+eta)];
        dNdxii  = .25*[      -(1-eta),       (1-eta),       (1+eta),      -(1+eta)];
        dNdetai = .25*[(1-xi)*-1     ,(1+xi)*-1     ,(1+xi)        ,(1-xi)        ];
        
        %% Jacobian
        Jac = ...
            [ [dNdxii;
            dNdetai]* obj.XE.';
            Ni*[obj.n1',obj.n2',obj.n3',obj.n4'].' ];
        
        detJ = det(Jac);
        
        %% Rotation matrix
        z_i = cross3(Jac(1,:),Jac(2,:))./norm(cross3(Jac(1,:),Jac(2,:)));
        x_i = cross3([0 1 0],z_i)./norm(cross3([0 1 0],z_i));
        y_i = cross3(z_i,x_i);
        T = [x_i;y_i;z_i];
        
        %% Partial derivatives wrt physical coordinates
        dNdxzy = T/Jac*[dNdxii; dNdetai; 0 0 0 0];
        %% Element thickness at integration point
        tg = t;
        
        % Calculate strain displacement matrix at xi & eta
        
        e1=[dNdxzy(1,1)     0           0;
            0               dNdxzy(2,1) 0;
            dNdxzy(2,1)     dNdxzy(1,1) 0];
        
        e2=[dNdxzy(1,2)     0           0;
            0               dNdxzy(2,2) 0;
            dNdxzy(2,2)     dNdxzy(1,2) 0];
        
        e3=[dNdxzy(1,3)     0           0;
            0               dNdxzy(2,3) 0;
            dNdxzy(2,3)     dNdxzy(1,3) 0];
        
        e4=[dNdxzy(1,4)     0           0;
            0               dNdxzy(2,4) 0;
            dNdxzy(2,4)     dNdxzy(1,4) 0];
        
        s1=[0 0 dNdxzy(2,1) 0    Ni(1) 0
            0 0 dNdxzy(1,1) Ni(1) 0    0];
        s2=[0 0 dNdxzy(2,2) 0    Ni(2) 0
            0 0 dNdxzy(1,2) Ni(2) 0    0];
        s3=[0 0 dNdxzy(2,3) 0    Ni(3) 0
            0 0 dNdxzy(1,3) Ni(3) 0    0];
        s4=[0 0 dNdxzy(2,4) 0    Ni(4) 0
            0 0 dNdxzy(1,4) Ni(4) 0    0];
        
        z3 = zeros(3);
         B  = [[e1  z3
        z3  e1
        s1                  ] * [T z3; z3 T*obj.A(:,:,1)],...
        [e2  z3
        z3  e2
        s2                  ] * [T z3; z3 T*obj.A(:,:,2)],...
        [e3  z3
        z3  e3
        s3                  ] * [T z3; z3 T*obj.A(:,:,3)],...
        [e4  z3
        z3  e4
        s4                  ] * [T z3; z3 T*obj.A(:,:,4)]];
    %% Modify constitutive (stress-strain) matrix for bending
    G = obj.G;
    ke =ke+ B([1:2,4:5],:).'*G([1:2,4:5],[1:2,4:5])*B([1:2,4:5],:) *tg*detJ;
    end
end
 % 1 point integration for shear
 xi=0;
 eta=0;
  %% shape function evaluations
        Ni      = .25*[(1-xi)*(1-eta),(1+xi)*(1-eta),(1+xi)*(1+eta),(1-xi)*(1+eta)];
        dNdxii  = .25*[      -(1-eta),       (1-eta),       (1+eta),      -(1+eta)];
        dNdetai = .25*[(1-xi)*-1     ,(1+xi)*-1     ,(1+xi)        ,(1-xi)        ];
        
        %% Jacobian
        Jac = ...
            [ [dNdxii;
            dNdetai]* obj.XE.';
            Ni*[obj.n1',obj.n2',obj.n3',obj.n4'].' ];
        
        detJ = det(Jac);
        
        %% Rotation matrix
        z_i = cross3(Jac(1,:),Jac(2,:))./norm(cross3(Jac(1,:),Jac(2,:)));
        x_i = cross3([0 1 0],z_i)./norm(cross3([0 1 0],z_i));
        y_i = cross3(z_i,x_i);
        T = [x_i;y_i;z_i];
        
        %% Partial derivatives wrt physical coordinates
        dNdxzy = T/Jac*[dNdxii; dNdetai; 0 0 0 0];
        %% Element thickness at integration point
        tg = t;
        
        % Calculate strain displacement matrix at xi & eta
        
        e1=[dNdxzy(1,1)     0           0;
            0               dNdxzy(2,1) 0;
            dNdxzy(2,1)     dNdxzy(1,1) 0];
        
        e2=[dNdxzy(1,2)     0           0;
            0               dNdxzy(2,2) 0;
            dNdxzy(2,2)     dNdxzy(1,2) 0];
        
        e3=[dNdxzy(1,3)     0           0;
            0               dNdxzy(2,3) 0;
            dNdxzy(2,3)     dNdxzy(1,3) 0];
        
        e4=[dNdxzy(1,4)     0           0;
            0               dNdxzy(2,4) 0;
            dNdxzy(2,4)     dNdxzy(1,4) 0];
        
        s1=[0 0 dNdxzy(2,1) 0    Ni(1) 0
            0 0 dNdxzy(1,1) Ni(1) 0    0];
        s2=[0 0 dNdxzy(2,2) 0    Ni(2) 0
            0 0 dNdxzy(1,2) Ni(2) 0    0];
        s3=[0 0 dNdxzy(2,3) 0    Ni(3) 0
            0 0 dNdxzy(1,3) Ni(3) 0    0];
        s4=[0 0 dNdxzy(2,4) 0    Ni(4) 0
            0 0 dNdxzy(1,4) Ni(4) 0    0];
        
        z3 = zeros(3);
         B  = [[e1  z3
        z3  e1
        s1                  ] * [T z3; z3 T*obj.A(:,:,1)],...
        [e2  z3
        z3  e2
        s2                  ] * [T z3; z3 T*obj.A(:,:,2)],...
        [e3  z3
        z3  e3
        s3                  ] * [T z3; z3 T*obj.A(:,:,3)],...
        [e4  z3
        z3  e4
        s4                  ] * [T z3; z3 T*obj.A(:,:,4)]];
    %% Modify constitutive (stress-strain) matrix for bending
    G = obj.G;
    ks = B([3,6:8]  ,:).'*G([3,6:8]  ,[3,6:8]  )*B([3,6:8]  ,:) *tg*detJ;
    ke = ke + 4*ks;
    %% transform to global coordinate system and save
    REG(22:24,22:24)=TEG;REG(19:21,19:21)=TEG;REG(16:18,16:18)=TEG;REG(13:15,13:15)=TEG;
    REG(10:12,10:12)=TEG;REG(7:9,7:9)=TEG;REG(4:6,4:6)=TEG;REG(1:3,1:3)=TEG;
    RGE = REG.';
    % save the I,J,K for assembly
    [ie,je,kke]=find(RGE*ke*REG);
    eid=k*ones(size(ie));
    I=[I;ELEMENT_DOF(k,ie)'];
    J=[J;ELEMENT_DOF(k,je)'];
    Kij=[Kij;kke];
    Eid=[Eid;eid];
end
clamped_nodes=find(COORD(:,3)==0);
clamped_DOF=[6*(clamped_nodes-1)+1;6*(clamped_nodes-1)+2;6*(clamped_nodes-1)+3;6*(clamped_nodes-1)+4;6*(clamped_nodes-1)+5;6*(clamped_nodes-1)+6];
load_nodes=find(COORD(:,3)==10);
extr_load_nodes=find(COORD(:,3)==10&(COORD(:,2)==10|COORD(:,1)==10));
extr_load_dof=[6*(extr_load_nodes-1)+3];
load_DOF=[6*(load_nodes-1)+3];
all_DOFs=(1:6*size(COORD,1))';
free_DOFs=setdiff(all_DOFs,clamped_DOF);
load_module=[ones(size(load_DOF));-0.5*[1;1]];
F=sparse([load_DOF;extr_load_dof],1,load_module,6*size(COORD,1),1);
K=sparse(I,J,E*Kij,6*size(COORD,1),6*size(COORD,1));
K=(K+K')/2;
U=zeros(size(F));
U(free_DOFs)=K(free_DOFs,free_DOFs)\F(free_DOFs);
figure
patch('Faces',ELEMENT,'Vertices',COORD,'FaceVertexCData',U(3:6:end),'FaceColor','interp');
colorbar
axis equal
view(-45,24)
%%
function p = cross3(u,v)
    p = [u(:,2).*v(:,3) u(:,3).*v(:,1) u(:,1).*v(:,2)]-[u(:,3).*v(:,2) u(:,1).*v(:,3) u(:,2).*v(:,1)];
end
##### SOURCE END #####
--></body></html>